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In this paper, we present an algorithm for the solution of the von Karman 
equations of elasticity theory and related problems. Our method of successive 
reconditioning is able to avoid convergence problems at any ratio of the non- 
linear streching and the pure bending energies. We illustrate the power of the 
method by numerical calculations of pinched or compressed plates subject to 
fixed boundaries. 
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The Navier-Stokes equation is not the only nonlinear partial differential equation in 
classical mechanics which has so far eluded a systematic analytical or numerical solution. 

Notoriously difficult partial differential equations appear also in the theory of elastic- 
ity; an example is given by the von Karman equations, which describe the elastic energies 
governing the deformation of a thin plate [[[]]. The von Karman equations are simpler than 
the Navier-Stokes equation in that they are local, though nonlinear. A variational theorem 
applies, which expresses the fact that the plate searches to minimize the total elastic energy. 

For a plate (described by the parametrization r(s,t) = (x,y, z)) of thickness h, the von 
Karman energy is given by [|TJ 
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where 
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Uij is the strain tensor with nonlinear terms in the deformations U\ = x(s,t) — s, u 2 = 
y(s,t) -t, z = z(s,t): 

1 dui duj 1 dz dz 

Oij is the stress tensor, linearly proportional to Uij. E e i and a are the Young modulus and 
the Poisson ratio, respectively. 

The appropriate framework for a precise numerical solution of this elastical problem is 
the finite element method 0, in which the plate (described by intrinsic variables s and t) 
is cut up into triangles, the corners of which come to lie at specified "nodes". Inside each 
triangle, any of the three functions is interpolated by a high-order polynomial, such that it 
is continuously differentiable as one passes from one triangle to the next one. This ensures a 
well-defined energy. Classic work on the finite element method has established the existence 
of a "magic" polynomial of order 5: 

x(s, t) = ax + a 2 s + a 3 t + a 4 s t + . . . + a 2 i£ 5 (2) 

(and similarly for y and z). The 21 parameters in eq. (Q) are fixed by prescribing at the 
three corners the function value x, the first and second derivatives dx/ds, dx/dt, d 2 x/ds 2 , 
d 2 x/dsdt, d 2 x/dt 2 as well as the normal derivatives dx/dn on the midpoints of the triangle's 
sides. We denote these variables of different origin by a symbol £ = (£1, . . . <^v), with iV the 
total number of variables. The polynomial eq. @ possesses an important symmetry property 
with respect to spatial rotations 0. It uses all the 21 polynomial parameters to satisfy, in an 
optimal choice, the continuity of the first derivatives across the sides of the triangle. Notice 
that all the second derivatives of the functions are imposed at the nodes. By construction, 
the functions x(s,t) etc, are thus twice continuously differentiable at these points. 

Given the nodes, and the variables £, it is straightforward to compute the interpolating 
polynomials (in the simplest fashion by inverting a 21 x 21 matrix), and to compute the 
local energy E(s(£), £(£)) as well as the total energy £(£) = / E(s(£),t(£))ds dt exactly 0. 



2 



As explained before, the finite element functions are valid variational test functions. It 
is therefore appropriate to search for the set of variables £ minimizing the total energy. This 
multidimensional minimization problem (solve for min^£{^)) is at the heart of our present 
concern, and our main result will consist in an algorithm which, for the first time, makes 
possible a direct attack at the von Karman equations in the presence of strong deformations, 
and at any value of the thickness |4]]. Any naive attempt to solve it is doomed to fail because 
of the large number of variables at hand. In addition, we have to solve the variational problem 
to great precision (essentially achieve |V£| = 0), since we are not really interested in the 
numerical value of the total energy, but in the geometrical aspects of the solution, which are 
much slower to converge than 8. 

A few algorithms are specifically geared at the solution of large minimization problems 
(for continuously different iable functions). They all attempt a local fit of the function by a 
parabola 

£{Z)~c + bt + ±ZHZ (3) 

where H is the Hessian matrix of second derivatives. Based on the knowledge of the function 
and of its numerical gradient b, strategies differ on how to economize on the computation of 

We have initially been extremely frustrated with the performance of these algorithms for 
large N, especially in the strongly nonlinear regime of small h. We illustrate the difficulties 
on a test example with iV = 272, a compressed half cylinder (c/fig. 2) of thickness h = 0.01, 
which will be further discussed later on. The upper curve in fig. 1 shows the total energy 
£ as a function of the iteration number, using one of the standard algorithms. For the 
given boundary condition, the simulation has proceeded for a few weeks on our work station 
without achieving reasonable convergence. In particular, the iteration never exhibits the 
quadratic convergence rate, which is the hallmark of the fast minimization algorithms, and 
which allows in principle the solution of problems with hundreds or thousands of variables 

Before exposing our solution for the present nonlinear case (i. e. when H depends on £), 
we shortly discuss the corresponding linear problem, in which the function H is independent 
of £. The behavior of iterative minimization routines, such as the Hestenes-Stiefel conjugate 
gradient method, has been discussed in the literature, (c/|| for a very clear discussion). The 
result is that the minimization algorithm performs well if the matrix H is well conditioned or 
is close to the unit matrix. Preconditioning algorithms have been devised, which transform 
the matrix H into a similar matrix, close to the identity . 

The fact that the nonlinear minimization program in fig. 1 initially works very badly 
can thus only mean one of two things: we either never enter the quadratic region, where 
the energy E{£) can be approximated by eq. (|3[), or we have to do with badly conditioned 
matrices H(£). It is an explicit computation of H (which is not normally undertaken), that 
has convinced us on several examples that the approximation eq. (|3|) becomes acceptable 
quite early in the simulation (it is to test this hypothesis that the very long initial runs 
were undertaken). However, the matrix H becomes very often extremely ill-conditioned 
especially if h is small. In the simulation followed in fig. 1, the eigenvalues of H span 6 
orders of magnitude at the point a) (specified by variables £ a ). 
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This has convinced us that a reconditioning of the matrix is necessary. We have used 
the simplest reconditioning possible: every so often, say at point a), we explicitly compute 
H (= H(£ a )) and its eigenvalues a, and eigenvectors rji. We then choose new variables 
£ — Vi/V®i- At point a), the matrix H thus ideally starts out as the unit matrix, and it 
then gets modified as we move away from £ a M. 

The expenditure of computing the matrix H(£ a ) and of reworking the complete calcu- 
lation in terms of £ may seem enormous. However, it is immediately rewarded by a strong 
decrease in the energy, and a concomitant fast motion in the variable space. This can be 
seen on the middle curve in fig. 1 starting at point a). As we move away from £ a , the con- 
dition of H (= H(£)) necessarily deteriorates again. In the example of fig. 1, starting from 
a well-conditioned matrix at point a), we reach at point b) an eigenvalue spectrum which 
again spans 5 orders of magnitude. In fact, most eigenvalues are of order 1, but a few very 
small eigenvalues corrupt the condition of the matrix. 
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FIG. 1. Evolution of the minimization in a test case (same problem as in fig. 2, h = 0.01, 16 
nodes, N = 272 variables). Upper curve: total energy £ vs iteration using a standard conjugate 
gradient algorithm. Reconditioning at a) and b) leads to the lower two curves. The insets show 
the local energy E(s,t) at a), b) and c). The final solution is given in fig. 2. 

If we continued our calculation, we would stay on the middle curve indicated in fig. 1, and 
undergo a gradual deterioration of the convergence. If, on the other hand, we recondition 
the matrix a second time, at point b), we immediately approach the solution of the problem 
(and witness a quadratic convergence rate (c/ ]5[]), which means that we double the number 
of significant places of the solution per iteration) . A zero gradient of £ (to machine precision) 
is reachable without problems. 

The final approach of the convergence is usually achieved after a few reconditionings, the 
precise number of which depends on the physical nature of the problem. A good criterion 
is to try reconditioning if £ no longer decreases even though |V£| is not yet approaching 
(which would indicate convergence). It is thus evident that reconditioning is useless in the 
initial stage of the minimization, in which even the standard algorithms decrease the energy 
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quite well. In fig. 1, we also show the local energy distribution, at points a), b) and c). 

The physical problem discussed so far consists in an elastic plate, whose zero-energy 
configuration is the unit square. We impose a cylindrical contour on two opposing sides. 
In addition, the plate is compressed along these two sides, as indicated in fig. 2. In our 
numerical work, we have to be concerned with the existence of local minima. In order to 
avoid the problem, we solve the minimization problem repeatedly for increasing values of 
the compression. For example, in the upper part of fig. 2, we first impose the cylindrical 
shape, and then compute minima of the von Karman energy for increasing values of the 
lateral compression. For all the values of the plate thickness ft, = 0.1, ... , 0.001 and of the 
compression, we have observed no numerical problems other than a gradual increase in the 
number of reconditionings required as h decreases. 




FIG. 2. Minimum energy solution for an elastic plate which was first bent, and then compressed 
(upper) and first compressed and then bent (lower). 

In the lower part of the fig. 2, we have followed a different procedure, by first compressing 
the flat plate, and then gradually imposing the cylindrical outer shape, in such a way that 
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the final constraints are exactly equivalent in the two parts of the figure. In both cases, we 
have undertaken extensive tests which have convinced us that precise numerical solutions 
are obtained with of the order of 16 independent nodes, and a few hundred independent 
variables ||. 

Surface Local Energy 




FIG. 3. Adaptive solution for the case of the pinched half-cylinder. The main figure presents 
contour plots of E(s, t) and the triangulation. 

Left two quadrants: Nonadaptive calculation with TV" = 271 variables, 16 nodes. 
Right upper quadrant: Nonadaptive, ./V = 126, 9 nodes 
Right lower quadrant: Adaptive, ./V = 126, 9 nodes. 

The adaptive calculation with 9 nodes gives essentially the same variational energy as the 
non-adaptive calculation with 16 nodes. Iterative minimization was used, as described in the 
text. 

Up to now, we have discussed the minimization problem only at fixed choice of the nodes. 
The most obvious way of increasing the precision of a calculation consists in adding new 
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ones, e.g. by subdividing existing triangles. This strategy has been explored in the literature 
|J, but soon encounters its limits as the number of variables increases very rapidly with 
the number of independent nodes. However, the node positions may be included in the 
variational problem. Let us denote the node variables by a symbol, 77. A sharper variational 
function can be obtained not only by increasing the dimension of rj, but by solving the 
extended minimization problem 

min V £E(£(r})) (4) 

Just as the Gaussian integration (in which the function to be integrated is evaluated at 
optimal positions), we can expect substantial gains in the quality of the variational functions. 
To show the usefulness of the approach, we present in fig. 3 the results of three separate 
numerical calculations on the problem of the half-cylinder which is pinched in the middle (we 
begin by compressing the elastic plate into a half cylinder, and then impose a fixed position 
for the center of the plate, cf. the upper left inset in fig. 3). On the left, we show a contour 
plot of the energy density, and the position of the nodes for a (non-adaptive) numerical 
calculation with N = 271 variables, and 16 nodes. On the right upper side, we show the 
corresponding results for a non- adaptive calculation with N = 126 variables, 9 nodes, and 
on the right lower side an adaptive calculation, with the same number of nodes. There are 
two important features of the adaptive solution: the global energy E is of the same quality 
as the much more expensive, non-adaptive solution with 16 nodes. Notice that the energy 
density resembles much more the one of the larger calculation, and that the nodes wander 
into regions of large local energy. 

In order to simplify the computation, we have, in fact not performed a full synchronous 
minimization over 77 and £, which seems unnecessary. In fact, it clearly appears that among 
the variables £, those which concern the function values, and the first derivatives converge 
much faster than the second derivatives (since the test functions are only once continuously 
differentiable). It is thus natural to suppose that the change of the node positions will 
have the largest influence on second derivatives. The minimization was done in an iterative 
form, in which the usual minimization (min^S(^) at fixed rj) is supposed to yield reasonable 
values for the functions x, y, z, and the first derivatives. At fixed function values and first 
derivatives, we then search for optimal values of rj and of the second derivatives, after which 
we again solve the usual problem, at the new rj. This completely rigorous procedure (every 
function encountered is a true test function of the von Karman energy) is then iterated 
several times. It quickly finds better nodes, which especially show smoother variations of 
the second derivatives of x, y, z across the triangle boundaries. 

In conclusion, we have discussed in this paper a very efficient method to solve numerically 
the von Karman equations. In the examples studied, our approach of successive recondi- 
tioning takes away all the convergence difficulties (in the different regimes of nonlinearity, 
i. e. h). Further work will have to show whether the reconditioning can be obtained at 
'reduced cost', without computing the full numerical Hessian. We also discussed the idea of 
optimizing with respect to the node variables. This adaptive choice of nodes was shown to 
be very useful. In the example, the nodes (and the sides of the triangle) migrate towards 
regions of very large local energy. We think it ultimately possible to take into account non- 
elastic terms in the von Karman equations, and thus produce a true numerical calculation of 
"crumpled paper" [|K|. It should be evident that the method may be of general usefulness 
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for non-linear minimization problems in high dimension. 
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